For each stimulus intenisty, what is the probability that the set of 12 repeated stimuli include zero?
To assess the probability, we modelled the data using the binomial test with a 50% probability of ‘success’ (positive rating arbitrarily chosen as success) at each stimulus intensity for each participant. Significance was assessed at the \(\alpha\) = 0.05 level, using a two-tailed probability.
# Import
data <- read_rds('data/SPARS_A.rds')
# Inspect
glimpse(data)
## Observations: 1,823
## Variables: 19
## $ PID <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID0...
## $ block <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A",...
## $ block_order <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1...
## $ intensity <dbl> 3.75, 1.50, 3.25, 1.50, 3.00, 2.75, 1.00, 2....
## $ intensity_char <chr> "3.75", "1.50", "3.25", "1.50", "3.00", "2.7...
## $ rating <dbl> -10, -40, -10, -25, -20, -25, -40, 2, -40, -...
## $ rating_positive <dbl> 40, 10, 40, 25, 30, 25, 10, 52, 10, 40, 54, ...
## $ EDA <dbl> 18315.239, 13904.177, 11543.449, 20542.834, ...
## $ age <dbl> 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, ...
## $ sex <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ panas_positive <dbl> 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, ...
## $ panas_negative <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
## $ dass42_depression <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ dass42_anxiety <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ dass42_stress <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ pcs_magnification <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,...
## $ pcs_rumination <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, ...
## $ pcs_helplessness <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
# Select columns
data %<>%
select(PID, intensity, rating)
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to sign
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x < 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'two.sided'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "SPARS A: Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | two-tailed p-value",
caption = "*'success' arbitrarily chosen as positive SPARS ratings",
x = 'Stimulus intensity (J)',
y = 'SPARS rating') +
scale_x_discrete(breaks = seq(1, 4, by = 0.5),
labels = sprintf('%.2f', seq(1, 4, by = 0.5))) +
scale_y_continuous(limits = c(-50, 50)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4) +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing = unit(0.1, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)
## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'SPARS') %>%
# Remove ratings of 0
# Ratings of 0 are uninformative
filter(rating != 0) %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
data %>%
group_by(PID, intensity) %>%
summarise(repeats = n())
## # A tibble: 63 x 3
## # Groups: PID [?]
## PID intensity repeats
## <chr> <dbl> <int>
## 1 FST1601a 2.25 12
## 2 FST1601a 2.5 12
## 3 FST1601a 2.75 12
## 4 FST1601a 3 12
## 5 FST1601a 3.25 8
## 6 FST1601a 3.5 11
## 7 FST1601a 3.75 9
## 8 FST1601a 4 12
## 9 FST1601a 4.25 12
## 10 FST1601b 2.25 11
## # ... with 53 more rows
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to sign
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x < 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'two.sided'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "SPARS B: Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | two-tailed p-value",
caption = "*'success' arbitrarily chosen as positive SPARS ratings",
x = 'Stimulus intensity (J)',
y = 'SPARS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(-50, 50)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)
## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'NRS') %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to whether
# the value is greater than 0 (minimum rating on NRS)
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x == 0,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'greater'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "NRS (zero: 0): Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | one-tailed (greater) p-value",
caption = "*'success' chosen as NRS rating values greater than 0",
x = 'Stimulus intensity (J)',
y = 'NRS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(0, 100)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
# Import
data <- read_rds('data/SPARS_B.rds')
# Inspect
glimpse(data)
## Observations: 2,149
## Variables: 10
## $ PID <chr> "FST1601a", "FST1601b", "FST1602", "FST1603", "...
## $ block <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ trial_number <int> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3,...
## $ intensity <dbl> 3.75, 2.75, 4.25, 3.00, 3.50, 2.50, 2.25, 4.00,...
## $ intensity_rank <int> 7, 3, 8, 3, 5, 4, 1, 8, 3, 3, 8, 3, 6, 4, 4, 8,...
## $ scale <chr> "SPARS", "SPARS", "NRS", "SRS", "SRS", "SPARS",...
## $ scale_combined <chr> "SPARS", "SPARS", "COMBINED", "COMBINED", "COMB...
## $ rating <dbl> 3, 1, 9, -78, -32, 17, 27, 4, 5, 7, 0, -69, 9, ...
## $ rating_eq <dbl> 3.0, 1.0, 4.5, -39.0, -16.0, 17.0, 13.5, 4.0, 5...
## $ rating_z <dbl> 0.5043367, 0.4190808, -0.3963035, -0.6689292, 0...
# Select data
data %<>%
# Extract trials rated using the SPARS
filter(scale == 'NRS') %>%
# Select columns
select(PID, intensity, rating)
# Check number of repeats at each stimulus intensity
#data %>%
# group_by(PID, intensity) %>%
# summarise(repeats = n())
# Nest data by PID and stimulus intensity
data_nest <- data %>%
group_by(PID, intensity) %>%
nest()
# Generate data
data_nest %<>%
# Add probability of success column
mutate(prob = 0.5) %>%
# Extract rating data from dataframe
mutate(data_vec = map(.x = data,
~ .$rating)) %>%
# Recode rating data as categories according to whether
# the value is greater than 0 (minimum rating on NRS)
mutate(data_cat = map(.x = data_vec,
~ ifelse(.x <= 15,
yes = 'negative',
no = 'positive'))) %>%
# Count the number of positive and negative ratings
## positive numbers arbitrarily listed first == 'success'
mutate(success_count = map(.x = data_cat,
~ c(length(.x[.x == 'positive']),
length(.x[.x == 'negative'])))) %>%
# Conduct binomial test (two-sided)
mutate(binomial_test = map2(.x = success_count,
.y = prob,
~ binom.test(x = .x,
p = .y,
alternative = 'greater'))) %>%
# Extract p-value from binomial_test
mutate(binomial_p.value = map(.x = binomial_test,
~ .x$p.value %>%
round(., 3))) %>%
# Categorise p-value using a p < 0.05 threshold
## Significant: distribution deviates significantly
## from the theoretical distribution
## No correction for multiple comparisons
## (too conservative for explorartory analysis)
mutate(significant_p.value = map(.x = binomial_p.value,
~ ifelse(.x < 0.05,
yes = 'yes',
no = 'no')))
For each paticipant, we plotted raw SPARS ratings at each stimulus intensity and colour-coded the data according to whether the p-value returned by the binomial test was significant (distribution of data points deviates significantly from the theoretical expected distribution).
data_nest %>%
# Select data columns
select(PID, intensity, significant_p.value) %>%
# Unnest data
unnest() %>%
# Join with original data
right_join(data) %>%
# Reclass intensity as an ordered factor
mutate(intensity = factor(intensity,
ordered = TRUE)) %>%
# Plot
ggplot(data = .) +
aes(x = intensity,
y = rating,
fill = significant_p.value) +
geom_hline(yintercept = 0) +
geom_point(colour = '#000000',
shape = 21,
size = 3,
alpha = 0.7) +
labs(title = "NRS (zero: 0 to 15): Binomial test of positive/negative rating distribution",
subtitle = "Probability of 'success' = 0.5* | alpha = 0.05 | one-tailed (greater) p-value",
caption = "*'success' chosen as NRS rating values greater than 15",
x = 'Stimulus intensity (J)',
y = 'NRS rating') +
scale_x_discrete(breaks = seq(1, 4.5, by = 0.25),
labels = sprintf('%.2f', seq(1, 4.5, by = 0.25))) +
scale_y_continuous(limits = c(0, 100)) +
scale_fill_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
scale_colour_viridis_d(name = 'Deviates significantly from expected distribution: ',
option = 'C') +
facet_wrap(~ PID,
ncol = 4,
scales = 'free_x') +
theme(legend.position = 'top',
panel.grid = element_blank(),
panel.spacing.x = unit(0.1, 'lines'),
panel.spacing.y = unit(0.5, 'lines'),
strip.text = element_text(margin = margin(t = 0.1,
b = 0.1,
r = 1,
l = 1,
'lines')),
axis.text.x = element_text(angle = -90,
hjust = 0,
vjust = 0.5))
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 magrittr_1.5 forcats_0.3.0 stringr_1.3.1
## [5] dplyr_0.7.6 purrr_0.2.5 readr_1.1.1 tidyr_0.8.1
## [9] tibble_1.4.2 ggplot2_3.0.0 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.18 cellranger_1.1.0 pillar_1.3.0
## [4] compiler_3.5.1 plyr_1.8.4 bindr_0.1.1
## [7] tools_3.5.1 digest_0.6.15 viridisLite_0.3.0
## [10] lubridate_1.7.4 jsonlite_1.5 evaluate_0.11
## [13] nlme_3.1-137 gtable_0.2.0 lattice_0.20-35
## [16] pkgconfig_2.0.1 rlang_0.2.1 cli_1.0.0
## [19] rstudioapi_0.7 yaml_2.2.0 haven_1.1.2
## [22] withr_2.1.2 xml2_1.2.0 httr_1.3.1
## [25] knitr_1.20 hms_0.4.2 rprojroot_1.3-2
## [28] grid_3.5.1 tidyselect_0.2.4 glue_1.3.0
## [31] R6_2.2.2 fansi_0.2.3 readxl_1.1.0
## [34] rmarkdown_1.10 modelr_0.1.2 backports_1.1.2
## [37] scales_0.5.0.9000 htmltools_0.3.6 rvest_0.3.2
## [40] assertthat_0.2.0 colorspace_1.3-2 labeling_0.3
## [43] utf8_1.1.4 stringi_1.2.4 lazyeval_0.2.1
## [46] munsell_0.5.0 broom_0.5.0 crayon_1.3.4